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Abstract. In this paper we present a new two-level iterative algorithm for 
tomographic image reconstruction. The algorithm uses a regularization technique, 
which we call edge-preserving Laplacian, that preserves sharp edges between objects 
while damping spurious oscillations in the areas where the reconstructed image 
is smooth. Our numerical simulations demonstrate that the proposed method 
outperforms total variation (TV) regularization and it is competitive with the combined 
TV -^2 penalty. Obtained reconstructed images show increased signal-to-noise ratio 
and visually appealing structural features. Computer implementation and parameter 
control of the proposed technique is straightforward, which increases the feasibility of 
it across many tomographic applications. In this paper, we applied our method to the 
under-sampled computed tomography (CT) projection data and also considered a case 
of reconstruction in emission tomography The MATLAB code is provided to support 
obtained results. 


1. Introduction 

Frequently, in X-ray computed tomography (CT), the amount of collected projection 
data is lower that it is required by the Nyquist sampling theorem [T] . In medical imaging, 
the restrictions are applied to minimize the ionizing radiation which can harm living 
tissue cells [2]. In material science, the aim is to better resolve temporal resolution via 
higher frame rate acquisition j3]. In such cases of limited data, iterative techniques can 
provide better reconstructions than analytical methods [1]. 

Dealing with ill-posed and ill-conditioned inverse problems, iterative methods 
require regularization to constrain the space of desirable solutions. Due to edge¬ 
preserving properties, total variation (TV) regularization [3] has been extensively used 
in tomographic iterative reconstruction for the last three decades. However, TV penalty 
produces piecewise-constant images (so-called “cartoon” or “staircasing” effect) even if 
the original object is smooth [3]. This “cartoon” effect can be particularly undesirable in 
emission tomography (ET) P], such as Positron Emission Tomography (PET) or Single- 
Photon Emission Computed Tomography (SPECT) where due to limited resolution 
of the imaging system and low photon-number statistics, reconstructed images are 
naturally blurred. Since in ET the activity distribution is piecewise-smooth and the use 
of TV penalty can result in undesirable artifacts, this can potentially bias the following 
clinical interpretation of reconstructed images [7]. 

Various improvements have been made to reduce “cartoon” appearance of the 
recovered images p n Eui [m [El na El. The most straightforward way is to couple 
TV term with £2 norm PEI [HE!. The goal of this approach is to establish a trade¬ 
off (by means of regularization constants) between piecewise-constant and piecewise- 
smooth regularization terms. This, however, complicates the optimization problem 
and for practical reconstruction purposes it can be a challenging task to establish the 
required regularization parameters. Therefore, the use of a single penalty term which can 
accommodate the required properties is highly desirable for tomographic reconstruction. 
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Some single term penalties have been proposed for image denoising and they based on 
the edge preserving Laplaeian [H] [Tl] or generalized forms of TV norm [IH] . 

In this paper, we propose a high-order penalty which has similarities with the 
Laplaeian based methods [m El bnt it also possess some of uniqne qnalities. We 
adapt our high-order penalty to the tomographic reconstruction problem using the two- 
step fixed point iteration algorithm ng. To demonstrate applicability of the proposed 
regularizer to image reconstruction problems we compare it with the classical TV [3] 
and combined TV-t '2 methods for CT and ET image reconstruction. All techniques are 
compared quantitatively and visually. 

2. Method 

2.1. Image reconstruction problem 

Tomographic image reconstruction problem consists of determining the inner structure 
of an object based on its X-ray observations from the several different angular positions. 
Incoming photons with different energies (due to absorption) are registered by detectors 
and information about the path length a photon has travelled along the line can be 
decoded [T]. Therefore, by solving the inverse problem where projection data is given, 
the degree of absorption or attenuation coefficient of the object can be recovered. In 
mathematical terms, the reconstruction problem can be formulated as the least squares 
(LS) problem: 

u = argmin ||Ati — 6 II 2 , (1) 

U 

where b e is a discrete function of the number of the detector bins and the 
observation angles describing the projection data (sinogram), u G is a function 
of spacial variables describing the observed object (e.g. density of the object’s material) 
and A : is a sparse system projection matrix (discrete approximation of 

the continuous Radon transform [T]) mapping the “space of objects” to the “space of 
observations”. In continuous space the operator A is an integral operator, and zero is a 
condensation point of its singular values, which makes the problem 0 ill-posed ng. 

Depending on the numbers of spacial grid points, detectors and angles, the matrix 
A can be “fat” (the number of rows is less than the number of columns) or “tall” (the 
number of rows is greater than the number of columns). Irrespective of the shape, the 
singular values of A form a very tight cluster near zero owing to the same property of 
the integral operator from which A derives [TH] . 

2.2. Regularization 

The quadratic functional J{u) = ||Atr — 6 II 2 with its “normal” form A*Au = A*b can 
be minimized by a suitable iterative minimization algorithm, e.g. Conjugate Gradient 
Least Squares (CGLS) algorithm [T7], however, the convergence of iterations can be very 
slow because of the poor conditioning of the Hessian A* A [13]. The slow convergence 
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of iterations is closely related to another difficulty in dealing with this kind of problem, 
which is the ill-posedness of the problem. Generally, the solution u of the problem Q is 
not unique (if A is “fat”), and even if it is, in practical computation u is indistinguishable 
from any w + if ||Afi ||2 is below the round-off error level, which, in the case at hand, 
may happen even if ||h ,||2 is substantial. 

Both difficulties can be tackled by a regularization technique, whereby Q is 
replaced with 

1 2 

Ma = argmin?/;„(w), = - ||Am - 6||2 aR{u), ( 2 ) 

u Z 

where R(u) is a suitable regularization functional, and a is a positive scalar parameter. 
The LS data term in 0 can be replaced by other fit-to-data functional, for 
example, Poisson log- likelihood functional, which is better suited noise model for ET 
reconstruction [3]. Optimization problem with regularized Poisson log-likelihood is 

Ua,p = argmax^„,p( 6 . Aw), 

U 


M 

pip, Aw) = ^{-[Aw]i -F bi log[Aw]i} -F ai?(w). (3) 

i 

Now we consider various regularization terms R{u) that can be used in (|^ and Q cost 
functions. A prime example of the regularization is known as Tikhonov’s [IS], where 

Re^{u) = \\u\\l. (4) 


Regularization makes the minimization problems Q and ([^ well-posed, the eigenvalues 
of the Hessian of ipap) being not less than a. Tikhonov regularization is quadratic, 
therefore high frequencies which are related to the object boundaries are penalized more, 
resulting in a smooth recovery of w. To preserve boundaries one needs to consider 
non-quadratic penalties, e.g. the TV penalty [Hj can significantly improve oscillatory 
solutions. The differentiable (due to small e constant) TV penalty is given as: 


Rpy(w) — |||Vw|e||^ — 




(5) 


Unlike Tikhonov’s penalty Q, TV can recover function while preserving sharp 
discontinuities, however, it also leads to the patchy appearance of the image [HI El 
M mi nzi Ha [H]. Since TV is unable to recover smooth variations alone, it can be 
coupled with ^ 2 -norm based penalty 0EIIIZIII31 resulting in the combined functional 


RTV-hi"^) — lll^'^l 


elll 


, ^ II I 
a 


2 ’ 


( 6 ) 


where fi is an additional regularization parameter for the quadratic term. In [Hj, the 
following combined functional has been proposed: 

(Auf 


RTV-i2p) — lll^'**|e|il + ~ 


|Vw| 


(7) 


where Aw = u^x + Uyy denotes an isotropic Laplacian. In this work, we will be using 
(0 and (|3 functionals for comparison with the proposed penalty. 
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It is intuitively obvious that the regularization parameter a in Q and (j^ must not be 
large. In order to get some further insight into the issue, let us estimate the difference 
between Ua and u assuming for simplicity that all singular values of A are positive. 

Let us assume R{u) = ||Rm|| 2 , where R is a square non-degenerate matrix, and 
denote M = A*A, and N = R*R. In the case at hand, Uq, (where a may be zero) 
satisfies 


(M + aN)w„ = A*6, 

which implies the following equation for the regularization error, = Ua — 
(M -h aN)ha = —aNti. 

Multiplication by yields 


((N + aNM-^N)ha, K) = -a(Nff, M-^Nha). 


-It 


Now, in the left-hand side of (10) we have 

((N + aNM-iN)ff„,ff«) > {'NK,h^) = (R*Rff„,ff«) = \\\RhX, 

and in the right-hand side (Nti, M“^Nh,Q,) = 

(NM-^Nff, h^) = (RM-^Nff, Rh^) < RM-^Nff \\Rh, 


ct||2 • 


Thus, (10) implies 

llRffalU < a 


RM-^Nm 


( 8 ) 

(9) 

( 10 ) 


( 11 ) 


2 . 4 . Edge-preserving piecewise-smooth regularization 


At first glance, the regularization appears to be merely a compromise move that distorts 
the problem so that it becomes solvable. While this is certainly so in the case of 
Tikhonov’s regularization, an alternative viewpoint can be offered, which is helpful 
in designing a proper regularization for the problem at hand. One observes that the 
problem 0 can only be solved approximately, if only because of the inexactness of 
computer arithmetic. One observes further that it may have infinitely many approximate 
solutions that are indistinguishable in practical computation, as pointed out in the 
previous section, if one is only guided by the smallness of the goal data fidelity functional 
||Ati — b|| 2 - The regularization can be viewed as some kind of additional criterion that 
helps to verify whether a particular computed solution is acceptable. This viewpoint is 
supported by the fact that in the case where R{u) = ||Rm|| 2 , the regularized problem 
(||) is equivalent to the original problem Q for these extended A and b: 


A 


a 


A 

aR 


b 

0 


( 12 ) 


In the problem (j^, A is a discretization of an integral operator, owing to which ||Ah ,||2 
is small on oscillating grid functions h with wavelengths that are close to the grid step. 
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Hence, if one directly applies e.g. CGLS algorithm to the minimization of the quadratic 
functional Q, then after sufficiently many iterations one is likely to end up with a 
quickly oscillating approximate solution ii. But most images that one encounters in 
practice do not feature such oscillations and can be actually represented by piecewise- 
smooth functions u. This suggests that the value of R{u) should be large on short- 
wavelength functions u. At the same time, R{u) should remain small on the boundaries 
(walls) between objects constituting the image, where u is discontinuous or has large 
gradients. 

The following regularizer is therefore suggested: 

2 2 

Rel{u) = 


d^u 

2 

d^u 


+ 

2 

dy^ 


where the weights Wx and Wy are given by 

-1 


w. 






-1 


(13) 


(14) 


/3 is a positive scalar parameter and it determines which points are considered to be on 
an edge between objects rather than inside it - the smaller the value of /3, the smaller 
the edge area. Constants and Oy are the average x- and ^-derivatives of u. These 
averages are introduced purely for the sake of scale-invariance, and can be computed 
e.g. as ttx = ‘^.Umax/dx and Oy = 2 umax/dy where Umax is the maximum of u and dx and 
dy are the sizes of the square containing the image. 

By design, Rel{u) (where EL stands for the Edge preserving Laplacian) is large 
on a short-wavelength functions u with small oscillation amplitude (approaching to 
isotropic smoothing) inside the reconstructed objects, and small on the edges between 
them (owing to the smallness of weights Wx and Wy). The edge points are identified as 
points where the derivatives of u are significantly higher than their average values. 

The regularizer ( [l^ shares some similarities with the high-order penalty proposed 
for image denoising in [TT| : 


u = argmin 

U 



A 

2 



u — uo)‘^dx, 


(15) 


where Uq is a noisy image and weights defined as uj = {1 -\-(3\VG„ * Uq\)~^ ^ where 
Ga denotes the Gaussian filter with the kernel width a and * is the convolution. 
There are few differences of the proposed penalty (13) from (15). For tomographic 
reconstruction a good initialization Uq is not usually available, therefore our regularizer 
is built on a different principle of “catching” prominent edges. The proposed penalty 
(13) is independent of the smoothed gradient while remain stable to the presence of 
noise (we will justify it with our numerical experiments). Moreover, we incorporate 
directional high-order smoothing in our term with consideration of variant weights 
(second derivatives are weighted independently), in our penalty Wx 7 ^ Wy (cross partial 
derivatives can be also added). Overall the penalty (13) is more rigorous than (15) and 
it will not allow any harmonic noise into the solution. 
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In order to simplify the computation of the gradient of and dealing with the 

minimization of a quadratic functional, we resort to an inner-outer iterative scheme 
with lagged diffusivity fixed point iteration [IS], with Wx and Wy only updated on 


restarts (this approach helps to deal with non-convexity of penalty (13) cf. widely used 
trust region technique dl)- 

Let us denote by La, and the matrices representing the discretized second partial 
X- and ^-derivatives with Neumann boundary conditions, and by and Wy the 
diagonal matrices representing and Wy. The symmetric gradient matrix R(ti), such 
as VR{u) = R{u)u, is given for EL term as (ignoring the dependence of Wx and Wy on 
u - cf. inner-outer iterative scheme) 

Rel = L-w;w.L, + qw;w„L,. (is) 

We compare the proposed EL penalty with TV regularization term ([^ with the following 
gradient matrix Rtv: 

Rtv = D:d>(tr)D, + D;d>(tr)D„ (17) 

where D^, and are matrices representing the discretized first partial x- and y- 
derivatives with Neumann boundary conditions, d>(ti) = diag{<fi{u)) is a diagonal matrix 
whose diagonal elements are fit) = 2\/t -\- e^. Different choices can be used for 

(/)(t) function to approximate £i norm, for example the Huber function [ISl El] • 


We also compare the proposed penalty (13) with TV -^2 functional with the gradient 
matrix defined as 


Rry-r, = D:T(tr)D, + D*T(tr)Dj, + KT{u)Lx + KT{u)l.y, 


(18) 


where \l/(w) is a diagonal matrix with elements Q;/|VM|e and T(m) is a diagonal matrix 
with elements 2/r/|Vw|^. Values e and 7 were chosen similar to 0 , 6 10 "^max 

'T = 

/ max * 

Here we used the lagged diffusivity fixed point iteration na (see algorithm to 
solve regularized LS optimization problem with penalties (IbpT) and 


Algorithm 1 Lagged Diffusivity Fixed Point Method for Regularized LS 

Uq := initialization 

begin outer (fixed point) iterations 

Qi, := A*{Au — b) -\- «Rtr; % gradient 
H := A*A -|- ckR; % approximate Hessian 
Sy := % quasi-Newton step 

:= Ui, -\- Si,; % update approximate solution 
n = n -Yl; 

end outer (fixed point) iterations 
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The system Hsj. = is solved by usual CG method ra. let I be an iteration 
index and L is the maximum number of inner iterations for CG method. Inner 
and outer iterations can be terminated earlier if the following stopping criteria met 
||s/ — Si+i ||2 < P and \\Ui, — Mi.+i ||2 < p respectively. We chose p = 10 to be a small 
constant. 

The largest eigenvalne of the matrix R is of the order 0{h~‘^), h = miii{hx,hy), 
and the smallest are of the order 0(1). If a is not very small (considerably larger 
than /j^), then the large condition nnmber of R*R can slow down the convergence of 
CG iterations for the minimization of 'Via(tt). To alleviate this problem, we introduce 
preconditioning that consists in the multiplication of the gradient of by the inverse 

of H = cr^I + aR in the course of CG iterations, where a is a scalar valne of the order 
of the largest singular value of A. Since matrix is very sparse, the application of its 
inverse to a vector can be efficiently performed (via the factorization of Ho-) by modern 
state-of-the art sparse direct solvers na. 

2.6. Regularized reconstruction with Poisson data fidelity term 

To solve regularized Poisson log-likelihood problem ([^ we use the splitting technique 
introduced in and used for SPECT reconstruction in 1211 . The main idea is to 
deconple data-fidelity term from the regnlarizer by solving two problems independently 
(see algorithm . From the structure of algorithm one can see that the optimization 


Algorithm 2 Splitting algorithm for regularized Poisson term 


tio := initialization with ones 
begin onter iterations 

i := u,./(A*l). * A*(b./Au,y, 

fo ■= 

I := 0 ; 


% MLEM step 


fi-\-i ■— fi + 'T{{fi — fo) + % image denoising 

:= /z+i; % update approximate solution 

u = u-\-l; 


end outer iterations 


problem with Poisson log-likelihood term is solved independently with Maximum 
Likelihood Expectation Maximization (MLEM) algorithm [IS]. Regularization is 
performed as an additional image denoising step. In MLEM step, 1 denotes the vector 
of I’s, and .* and ./ denote componentwise multiplication and division, respectively. 
Additionally the weights Wx and Wy for the proposed term (16) are calculated after 
each MLEM step and kept fixed in image denoising npdates. 
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In this section we present two different numerical experiments with the proposed 
penalty EL (16) and compare it with TV ( |I7| ) and TV-t '2 (18) penalties. In the hrst 
experiment we model reconstruction of the synthetic phantom with algorithm [T] and in 
the second experiment we perform reconstruction of more realistic phantom for emission 
tomography with algorithm The MATLAB code for tomographic reconstruction using 
TV ( [Tfl ) regularization and the proposed EL penalty (16) is provided ng. 


3 .1. CT reconstruction 

To test the proposed penalty in regularized tomographic reconstruction we designed an 
analytical phantom which consists of a smooth (two Gaussians and two parabolas) and 
piecewise-constant (one rectangle) functions (see Eig. [^. 



Figure 1. (a) 2D analytical phantom which consists of 2 parabolas, 2 Gaussians and 
a rectangle; (b) a surface representation of the phantom. 


To avoid of reconstructing on the same grid where projection data was generated 
(so-called reconstruction with “inverse crime ” [E]), we used a higher resolution of the 
phantom on a 500 x 500 isotropic pixel grid to generate projections with a strip kernel 
[T]. Then Poisson distributed noise was applied to the projection data, assuming an 
incoming beam intensity of 3T0^ (photon count). Reconstructions were calculated on 
a 250 X 250 isotropic pixel grid and with a linear projection model [I]. We used 90 
projection angles in 180 degrees (assuming a parallel beam geometry) to reconstruct the 
phantom. 

For a fair comparison of different regularizing penalties we initially optimized the 
regularization parameters (see Fig. with respect to the value of root-mean-square- 
error (RMSE), defined as: 


I'U 

— 

u\ 

I 2 


\u\ 

12 


(19) 


where u is an exact image and ti is a reconstructed image. 


We found empirically that /3 = 0.03 for EL penalty (16) gives good results for the 
presented experiments, therefore we will keep it fixed for the rest of our tests. With 
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Optimization for CGLS-TV reguiarization parameter 



a 


Optimization for CGLS-T\/L2 reguiarization parameter 



u 


Optimization for CGLS-EL reguiarization parameter 



a 


Figure 2. Optimization procedure to find the optimal regularization parameters for 
TV (left), TV -^2 (middle) and EL (right) algorithms. For TV -^2 penalty we optimized 
the second regularization constant (/i) while keeping the optimal a fixed (found for 
TV previously). 


fixed optimal regularization parameters (see Fig. we perform L = 80 outer (fixed 
point) iterations and 5 inner iterations of algorithm with different penalties (see Fig. 
[^. In Fig. one can see that the CGLS method diverges quickly while regularized 


Convergence plots for CGLS, CGLS-TV, CGLS-TVL2 and CGLS-EL methods 



Figure 3. Plot of RMSE values to the number of iterations for CGLS, CGLS-TV, 
CGLS-TV-f ?2 and GGLS-EL methods. GGLS-EL method outperforms other algorithms 
(see table [^. 

methods converge to a steady-state solution. The lowest RMSE value is achieved with 
the proposed EL penalty (see table [^. Reconstruction with TV penalty gives the 
highest RMSE value (from all compared regularized methods) however it shows faster 
convergence on first fixed point iterations. TV -^2 penalty performs slightly better than 
TV, but still with higher RMSE than EL method. 


Table 1. RMSE values for GGLS, GGLS-TV, GGLS-TV-tg and CGLS-EL methods 



CGLS 

CGLS-TV 

CGLS-TV-G 

CGLS-EL 

RMSE 

0.1710 

0.0919 

0.0902 

0.0868 


Reconstructed images are presented in Fig. Since CGLS-TV -^2 reconstruction 
might look more appealing than CGLS-EL we also show the surface representations of 
reconstructed images (see Fig. and horizontal middle cross-sections (see Fig. |^. 
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Figure 4. Reconstructed phantoms with (a) CGLS, (b) CGLS-TV, (c) CGLS-TV -^2 
and (d) GGLS-EL method. 



Figure 5. Surface representations of the reconstructed phantoms with (a) GGLS, (b) 
GGLS-TV, (c) GGLS-TV -^2 and (d) GGLS-EL method. 



Figure 6. Horizontal middle cross-sections of the phantom and reconstructed images 
(see Fig. |^. Left side images show the region selected and the right side images show 
the magnified region. It can be clearly seen that the EL penalty performs much better 
for smooth objects while slightly more perturbed in uniform areas (e.g. the top of the 
rectangle). 


One can notice that CGLS reconstruction is very noisy. CGLS-TV method better 
suppresses noise, however smooth features are strongly affected by the “staircasing” 
effect. CGLS-TV -^2 method provides reconstruction with smoother features and GGLS- 
EL method resolves smooth features even better (e.g. cone-shaped parabola). Although 
CGLS-EL method performs very well for smooth objects one can notice the wave-like 
variations of intensity in the background and also at the top of the rectangle (see Fig. 




















































Direct high-order edge-preserving regularization for CT 


12 


[^. This issue can be explained by the properties of our regularizer, in contrast to TV, 
our penalty does not seek the sparsest solution and does not penalize strongly (pushing 
to the constant value) a small intensity perturbations. The EL term tends to preserve 
all sharp edges while uniform noise is smoothed isotropically with the Laplacian. In Fig. 

one can see that the CGLS-EL method provides better recovery of smooth features 
while slightly higher (compare to TV and TV-t' 2 ) perturbations visible in uniform areas 
(the top of the rectangle), however, the edges of the rectangle are defined sharper with 
the EL penalty. 

3.2. ET reconstruction 

To simulate emission tomography reconstruction we designed a more realistic phantom 
from the high-quality X-ray scan of a mice bone. The data was acquired on a Nikon 
Metris Custom Bay cone-beam scanner at the Henry Moseley Manchester X-ray facility, 
and was reconstructed with the Feldkamp algorithm (see Fig. 0(1 eft)). We thresholded 
the obtained reconstruction and added six gaussians with various kernel widths (see Fig. 
(middle and right)). 



Figure 7. (left) High quality X-ray reconstruction of a mice bone; (middle) phantom 
created from the reconstruction (left) with a six Gaussians added to it; (right) a surface 
representation of the phantom. 


To simulate PET projection data we used NiftyRec [22], a software for tomographic 
reconstruction, providing GPU-accelerated reconstruction for emission and transmission 
computed tomography. The phantom size is 400 x 400 pixels and 300 projections was 
simulated. Poisson noise was added to projections with an expected number of 10 • 10® 
photon counts in total. Twenty noise realizations were simulated to estimate methods 
quantitatively. The point spread function of the PET system was modelled (with 
convolution of the sinogram columns with a Gaussian of full width half maximum of 
three pixels) in the projection and back-projection operations. No scatter was simulated 
in this study. For our experiments (see algorithm we performed 130 MLEM iterations 
and 5 inner iterations (denoising step). 

To quantify obtained reconstructions we used averaged over all noise realizations 
RMSE (19) values in the bone region (BR) and in Gaussian regions (GR). All 
regularization parameters were carefully selected by comparing the mean of all RMSE 
values over all noise realizations in GR and BR (see Fig. |^. 
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cc 

OQ 

LU 

LO 

z 

Qc: 


Optimization for TV reguiarization parameter (ET case) 


Optimization forTV-L2 reguiarization parameter (ET case) 


Optimization for EL reguiarization parameter (ET case) 




Figure 8. Optimization procedure to find the optimal regularization parameters 
for MLEM-TV (left), MLEM-TV -^2 (middle) and MLEM-EL (right) algorithms. Eor 
MLEM-TV -^2 method we optimized the second regularization constant (/i) while kept 
the optimal a fixed (found for MLEM-TV method previously). 


After estimation of regularization parameters we performed twenty reconstructions 
for each method with various Poisson noise distributions. The mean values for GR and 
BR over all noise realizations are shown in Fig. This result proves that the EL 
penalty is very successful in resolving smooth features (six Gaussians in this case) and 
also quite competitive for the BR (lower RMSE value than for TV). 


CC 

OQ 

LU 

cn 


Best RMSE values for MLEM, EM-TV, EM-TVL2 and EM-EL methods 

0.06 
0.059 
0.058 
0.057 
0.056 
0.055 
0.054 
0.053 
0.052 
0.051 

0.049 0.05 0.051 0.052 0.053 0.054 0.055 0.056 0.057 

RMSE GR 


Figure 9. Plot of the best RMSE values for MLEM, MLEM-TV, MLEM-TV-^g and 
MLEM-EL methods. MLEM-EL method outperforms other algorithms for GR and 
gives better result than TV for BR. 


we demonstrate reconstructed images for the best selected RMSE 

values. 


In Fig. [TWand 11 



Figure 10. Reconstructed phantoms with (a) MLEM, (b) MLEM-TV, (c) MLEM- 
TV -^2 and (d) MLEM-EL method. 


In Fig. and one can notice that the BR is very smooth for TV and TV- 
£2 penalties and some long-wave oscillations can be seen in the reconstructed image 
with EL penalty. This result corresponds to the expected behaviour of the EL penalty. 
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Figure 11. Surface representations of the reconstructed phantoms with (a) MLEM, 
(b) MLEM-TV, (c) MLEM-TV -^2 and (d) MLEM-EL method. 


We note here that the phantoms backgronnd (see Fig. is not as flat as TV and 
TV-t '2 penalty recovered it. Fnrthermore, a small size dot-like featnre (approximately 
in the centre of the phantom) is almost smoothed out with TV and TV -£2 recovery. 
However, it is visible and well recovered with EL penalty. The sharp features, overall, 
are reconstructed very well with MLEM-EL method and seem even sharper compare to 
other methods (see the bone outer rim in Eig[Iol). 

4. Discussion 

In this paper, we emphasize that piecewise-smooth images are more favourable than 
piecewise-constant ones, this, however might not be the case for all reconstructed 
objects. Therefore, some prior knowledge about the investigated object is needed to 
choose an appropriate regularization term. For instance, the activity distribution in ET 
is smooth, therefore, penalties like EL are particularly suitable. The proposed penalty 
gives more gentle approach to regularization in uniform areas and does not penalize small 
perturbations aggressively. Moreover, the proposed method can recover small features 
successfully (e.g. lesions in ET case) with both step or ramp intensity variations. 

In terms of the choice of parameters, complexity of computer implementation and 
the speed of computation, our method is very similar to reconstruction techniques with 
TV penalty. Our future work will be to explore further the space of parameters and 
give some recommendations for automated choice of some parameters. Additionally we 
will look into the issue of low-frequency oscillations of our penalty in the uniform areas. 

5. Conclusion 

In this paper, we presented a novel two-step iterative reconstruction algorithm with 
high-order regularization penalty. Our method outperforms in terms of signal-to- 
noise ratio the conventional total variation based reconstruction techniques and it 
is competitive with other state-of-the-art high-order based approaches. From the 
preliminary experiments, the proposed method is well suited for the limited data 
problems in X-ray tomography as well as emission tomography. 
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